CORRELATIONS FOR THE NOVAK PROCESS 



ERIC NORDENSTAM AND BENJAMIN YOUNG 



Abstract. We study random lozenge tilings of a certain shape in the plane called the 
Novak half-hexagon, and compute the correlation functions for this process. This model was 
introduced by Nordenstam and Young (2011) and has many intriguing similarities with a 
more well-studied model, domino tilings of the Aztec diamond. The most difficult step in 
the present paper is to compute the inverse of the matrix whose (i, j)-entry is the binomial 
coefficient C(A, Bj — i) for indeterminate variables A and 2?i, . . . , B n . 

Nous etudions des pavages aleatoires d'une region dans le plan par des losanges qui 
s'appelle le demi-hexagone de Novak et nous calculons les correlations de ce processus. Ce 
modcle a ete introduit par Nordenstam et Young (2011) et a plusieurs similarites des pavages 
aleatoires d'un diamant azteque par des dominos. La partie la plus difficile de cet article est 
le calcul de l'inverse d'une matrice ou l'element est le coefficient binomial C{Bj — i, A) 
pour des variables A et B\, . . . , B n indetermines. 

1. Introduction 

This paper is a continuation of the work in |NY11] . in which we initiated a study of 
the Novak half-hexagon of order n. This is a roughly trapezoid-shaped planar region (see 
Figured]), which can be tiled with 3n(n+l)/2 lozenges — rhombi composed of two equilateral 
triangles. The number of these tilings is computed in |NYllj to be 2 n<n+1 )/ 2 , the same as 
the well-studied Aztec diamond (see |EKLP92] ) and possesses a domino shuffling algorithm 
closely related to that of the Aztec diamond. We were able to exploit this similarity to 
prove an "arctic parabola" -type theorem for the Novak half- hexagon: that with probability 
tending to 1 as n — > oo, the tiling is trivial exterior to a parabola tangent to all three sides 
of the figure. 

The power-of-two tiling count, the existence of a domino shuffle and the simple limiting 
shape strongly suggest that it will be tractable to carry out the usual "next step" in the 
study of random tilings: namely, computing correlation functions for the tiling. Loosely 
speaking, the fc-point correlation function gives the probability that a fixed set of k lozenges 
will all be present in a lozenge tiling chosen with respect to the uniform measure on the set 
of all 2 n ( n+1 )/ 2 such tilings. There are a number of ways to compute these probabilities, all 
of which rely on the fact that the correlation functions are determinantal, meaning that they 
can be computed as the determinant of a k x k matrix, whose entries are evaluations of a 
correlation kernel. 

If these probabilities can be computed exactly, one can attempt to do asymptotic analysis 
of the correlation functions, and demonstrate that the tiling exhibits universal behaviour. 
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Here, universal is a loaded, technical term coming from statistical mechanics and random 
matrix theory: it means that the correlation functions tend to one of a handful of well-studied 
and frequently-occurring limit laws which originally come from random matrix theory. For 
instance, at points near the "arctic parabola", the correlations should tend to the Airy 
kernel (see [Joh05aJ) and in the bulk, they should tend to the Sine kernel. Many point 
processes exhibit these limit laws and other related ones, including eigenvalue distributions 
of random matrices [For 10] . the Schur process |UR03j . the length of the longest row of a 
random permutation OkoOO. B DJ99] . continuous Gelfand-Tsetlin patterns [Met 11] . domino 
tilings of the Aztec Diamond [Joh05aJ, lozenge tilings of the regular hexagon [Joh05b] and 
many more. 

1.1. Results. In this paper, we compute the correlation kernel for a rather general class of 
lozenge tiling problems, of which the half-hexagon is one (we cannot say anything about its 
asymptotics yet). The starting point of our method is the Eynard-Mehta theorem, explained 
in Section [3J This is a rather general theorem for computing the correlation functions for 
processes which can be described as a product of row-to-row transfer matrices, as ours 
can. The Eynard-Mehta theorem gives the correlation kernel in terms of the inverse of a 
certain matrix M. For the half-hexagon, M turns out to be the Lindstrdm-Gessel-Viennot 
matrix [Lm73l IGV85] . 

n + 1 

which computes the number of tilings of the order- n half-hexagon. In fact, our methods 
required us to invert a much more general matrix. 

Theorem 1. If A, Bi(l < i < n) are parameters and 
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Then, the Eynard-Mehta theorem yields the following corollary, which will be shown in 
Section [3J 

Corollary 2. The correlation functions for the Novak half-hexagon are determinental, with 
kernel given by 
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Figure 1. To the left is a random tiling of an order 20 Novak half- hexagon. 
To the right is the same tiling but rotated and skewed. As shown, a tiling 
corresponds to a set of non-intersecting random walks that at each time-step 
either stay or jump one unit step. 



1.2. Inverting a matrix. Inverting a fixed matrix of numbers is trivial in a computer. 
Symbolically inverting an infinite family of matrices with many parameters is much harder 
and comprises the bulk of the work in this paper. 

We inverted M with Cramer's rule: compute the adjugate matrix Aji (the transposed 
matrix of cofactors) and divide by the determinant of M. Krattenthaler [Kra99] gives many 
methods of evaluating such determinants; indeed, his Equation (3.12) allows us to compute 
det M. Computing the determinant of the matrix of the adjugate matrix, however, is signifi- 
cantly harder, so we first guessed the answer using the computer algebra system Sage [S + llj . 
The manner in which this guessing was done was itself nontrivial and may be of interest to 
others trying to invert matrices; some details are given in Section [21 

Once we had conjectured the form of Theorem [1] and simplified it considerably, we were 
able to prove it simply by showing that MM -1 is the identity matrix. 

1.3. Related Work. Metcalfe |Metllj has developed an alternative approach to problems 
of this type, by developing a theory of the asymptotics of a sort of interlacing particle process. 
The theory currently covers a slightly different setting, in which the positions of the particles 
is continuous, but Metcalfe is in the process of extending his methods to the discrete setting. 

A natural extension of this procedure would be to apply the ideas of Borodin- Ferrari [BF08J 
to analyse the dynamics of the domino shuffling algorithm described in |NYllj . 

In |Joh05b| . there appears a slightly less general kernel, written in terms of the Hahn 
polynomials; this is used to prove some theorems on the fluctuations of the frozen boundary 
of lozenge tilings of a hexagon. 
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2. An inverse matrix 



Recall that we want to compute the inverse of the matrix M from ([2]) by computing co- 
factors. The method of computation is the standard approach of experimental mathematics: 
First we guess the answer, making no attempt to be mathematically rigorous. Then, we 
prove our guess rigorously, by showing that MM -1 is the identity matrix. As the reader 
may imagine, the proof alone is not too helpful for guiding people who want to tackle similar 
matrix inversions in their own work, so we include here an account of how we were able to 
guess the expression in Theorem HJ 

Since M is symmetric in its columns, striking out column k simply removes all instances 
of the variable B^ from M. Rename the remaining B- variables B\, . . . , B n -\. Removing a 
row, say number s, is more complicated and splits the matrix into two blocks. To get ready 
to make our first round of guesses, we take out as many factors as possible so that the matrix 
elements are now integer polynomials in the variables. The remaining matrix can be written 
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Let P n , s (A, B) be the value of the second determinant. Because P nyS is antisymmetric in Bi, 
it is divisible by the order n — 1 Vandermonde determinant; once this is done, the remaining 
portion is symmetric, so we expand it as a (linear!) combination of the elementary symmetric 
functions e\. We started by computing P for a few different values of the parameters n and 
s. For s — 1 one quickly conjectures 
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where A means taking the Vandermonde determinant in the variables. For s — 2, sage gave 
us 



P 3>2 (A, B) =(A + i)A(B)(—2e 2 (B) + (A + A)e 1 {B) - (3A + 8)) 

P 4)2 (A, B) ={A + 1)\A + 2)A(B)(-3e 3 (B) + (A + 6)e 2 (B) - (3A + 12)ei(B) + (7 A + 24)) 
P 5>2 {A, B) ={A + If {A + 2) 2 {A + 3)A(B)(-4e 4 (S) + {A + 8)e 3 {B) - (3 A + W)e 2 {B) 

+ (7 A + 32)ei(S) - (15 A + 64)) 
P 6t2 (A, B) =(A + l)\A + 2f(A + 3)\A + 4)A(B)(-5e 5 (B) + (A + 10)e 4 (S) 

- (3A + 20)e 3 (5) + (7 A + 40)e 2 (5) - (15 A + 80)e 1 (B) + (31A + 160)) 



Following the immortal advice of David P. Robbing, we wrote the coefficients of this four- 
parameter expression in a tidy fashion, and applied the standard tools in experimental 
mathematics [OElTTl IWikllj to all the integer sequences we noticed. There were many 
patterns. For instance, the Stirling numbers of the second kind S(n, k) appeared in some 
the coefficients, as did the numbers n k and (n + l) k — n k . Since the Stirling numbers have 
the form 



k 



3=0 



and since all of the coefficients we computed seemed to grow exponentially as the index of 
the elementary symmetric function / decreased, we made the following ansatz: 

Ansatz 1. The coefficient of A k e n ^i_i(B) in P n s is of the form 

1 s 

-7 f*--i -J 

S - 3=0 

where fk,i,s,j is a low-degree polynomial. 

We asked sage to find polynomials fk,i,sj i n Ansatz [TJto fit the data, and to compute more 
terms. Computing more terms required heavy optimization of the sage code and, eventually, 
running the code on a very powerful computer. After once again writing fk,i, s ,j{n) in a tidy 
table and dividing out some obvious common factors, we noticed a new set of patterns: some 
of the fk,i, s ,j(n) were zth derivatives of the falling factorial functions (n — l)(n — 2) ■ ■ • (n — k). 
As such, we made a second ansatz: 

Ansatz 2. All of the fk,i,s,j(n) o,re linear combinations of falling factorials or their deriva- 
tives. 

Again, we asked Sage to compute the coefficients of these linear combinations for the data 
we had. This time we were able to guess the formula completely. In the end we conjectured 
that 



n-2 

(8) P n>s (A,B) = A(iv)H(A + r) n - 1 - r x 

r=l 



5-1 S3 Ak -l ( fb\ / / i \ i 

xEEEB-D^^f.(.-H,H (!) (»-u ■••<»- 

1=0 k=0 3=1 i=0 v ' \ v 7 

where s(n, k) are the Stirling numbers of the first kind. 

Obviously, ([5} needs to be simplified. By the generating function for the Stirling numbers, 
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When faced with combinatorial enumeration problems, I have a habit of trying to make the data look 
similar to Pascal's triangle" . |Rob91] 



By the Binomial Theorem, 
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by the definition of Binomial coefficients, 
(11) {n + A-l)---{n + A- ])= j\ 
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and lastly, by the generating function for the elementary symmetric polynomials 
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With these simplifications we can write 
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Now to get the inverse matrix we should transpose the cofactor matrix and divide with 
the determinant of the full matrix. The latter can be found through 
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which is a special case of [Kra99l equation (3.12)]. After a bit of simplification, Cramer's 
rule then leads us to conjecture that ([3]) is the inverse of M. 

of TheoremUi We have now, through these computer experiments, found a formula which 
we believe expresses M~ l . To prove that this guess is correct, we need to show that either 
MM -1 — I oi that M _1 M = / using that formula. One of these (the latter) is easy, the 
other is hard. First we write 



(15) [MM- 1 ]^ = Y,[M}aA M ~ V 



Next, we need the following technical lemma, to remove the variables B{ from the equation. 
Lemma 3. 

(16) £U, s -iJ U-JJL^rl *-l J U-7 

Proof. Recall from an undergraduate course how Lagrange interpolation works. Let's say 
you want to fit a polynomial y = p(x) of degree n — 1 to points (xi, yi), . . . , (x n , y n ). What 
you do is you define functions 



(17) Mx)= n 
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and then you compute your polynomial p by 

n 

(18) p(x)=^2y P Tp(x). 

The sum in the LHS of the Lemma is of exactly this form. Moreover, 

-(t - 1) • • • (t - a+ l)(A - t + n) ■ ■ ■ (A - t + a + 1) 
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is a polynomial of degree a — 1 + n — a = n — lint. So this sum does Lagrange interpolation 
of degree n — 1 to an expression that is already a polynomial of that degree. Replacing the 
sum with the correct polynomial proves the Lemma. □ 

Application of Lemma [3] reduces (IT5|) to 



This sum can be computed through Vandermonde convolution, as in [GKP891 Equation 
(5.25)], showing that 

[MM _1 ]q, i7 = (°_J = Sa n , 

which proves that we have indeed found the correct inverse matrix. □ 

3. Eynard-Mehta theorem 

In order to compute correlation functions, we must first describe tilings of the Novak 
half- hexagon as an ensemble of nonintersecting lattice paths (see Figured]). 

Consider n walkers on the integer line, started at time at positions x\ , x% , • • ■ , 
At time iV they end up at positions x 2 , • • • , At tick t of the clock they each 

take a step according to the transition kernel 4> t - In our special case, they either stay where 
they are or move one step to the right: 

(19) (f) t (x,y) = 5 X)V + 8 x+ltV , t = 0, . .. , N - 1. 

In addition, they are conditioned never to intersect. Let the positions of the walkers at 
time t be denoted = (x® < ■ ■ ■ < in') G N n and let a full configuration be denoted 
x = (x(°),...,xW). 

Then uniform probability on these configurations can be written 

N-l 

(20) p(x) = -l[dBb[Mx?\xf ¥1 %^. 

t=0 

The normalization constant Z is the total number of configurations. For the sake of 
notation define the convolution product * by 

/ * g(x, z) = ^2 f(x, y)g(y, z) 
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and let 

*-~*<fH-i)(x,y), s<t 



Vrs 

0, otherwise. 



By the Lindstrom-Gessel-Viennot Theorem [Lin73[ ICTV 85J. the total number of configura- 
tions is given by the determinant of the matrix 

(21) M=[<p , N (xf\xf ) )]Z =l . 

Correlations can now be computed using the Eynard-Mehta Theorem. Readable intro- 
ductions to it can be found in [ForlO] Section 5.9] as well as in [BR05J. 

Theorem 4 (Eynard-Mehta). Let m be a positive integer and let (t\,xi), . . . , (t m ,x m ) be a 
sequence of times and positions. The probability that there is a walker at time U at position 
Xi for each i = 1, . . . , m is given by 

det[A'(t i ,x i ;t i ,x i )]S=i 
where the function K , called the kernel of the process, is given by 

n 

K(r,x;s,y) = -<f> r ,,(.x,y) + ^ (j) r)N {x,xf f) )[M''\ j (t) QiS {x) °\y) 

In our particular case the walkers are going to start densely packed. At first we shall leave 
the end time N and the endpoints unspecified, i.e. 

(0) 

X\ = I, 
(iV) 

x\ = Vi, 

for i — 1, . . . , n. The particular transition function (fT9"|) gives </v jS as defined in (J5J). Inserting 
that into (|2ip gives 

N 



M 



which is exactly the matrix we inverted in the previous section. The kernel can then be 
written 

(22) K(r, x; s, y) = -(j) r , s (x, y) 

, y & (A) y (N + n - 1\ (N - 1 + j - fc\ k+1 n k-m_ 

h ( N :x) h\ k ~ i a 3-k r ] Jhvi-Vi- 

We state the result in this generality because the kernel derived in |Joh05bj is a special case 
for suitable choices of iV and yi in the sense that they are correlation kernels for the same 
process. It is not at all clear how to algebraically relate ( 122]) with the formula in |Joh05b[ 
Theorem 3.1], since the latter is a sum involving products of Hahn polynomials. 

In our particular case N = n + 1, and the end positions are fixed as y^ = 2i for % — 1, . . . , 
n. This specialisation leads to the expression in Corollary |2J 
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